// bandprocessing.h
// Library of spectral processing functions. 
// Author: Garreth Ruane

#include<stdio.h>
#include<math.h>

/* ------------------------------------------------------------------------------
List of spectral processing functions:

 calcAndDisplayBandParameters - Given a spectrum and a list of continuum boundaries, the band parameters are calculated and printed to the command line for each band. 

 calcParameters - Saves band minimum, band depth, FWHM, and Asymmetry to pointers passed to calcParameters function. 
 
 calcAsymmetry - returns asymmetry of band.
 
 removeContinuum - removes continuum, saves the continuum removed data to newdata.
 
 calcFWHM - Returns the Full Width at Half Maximum of band. 
 
 calcFWHMpts - Saves two coordinates interpolated to fall on the half max value. 
 
 calcBandDepth - Returns the difference between continuum and band minimum.
 
 calcBandMin - Returns wavelength at the band minimum based on parabolic fit.  
 
 
 
List of utility functions:

 parabolaEval - Evaluates a parabola function of the form c2x^2 + c1x + c0 @ xval.
 
 calcparabolaMinx - Returns x-coordinate of minimum of parabolic function.
 
 determinant3by3 - Returns the determinant of 3x3 matrix.
 
 copy3by3 - Copies 3x3 matrix "in" to 3x3 matrix "out".
 
 solve3by3system - Assigns the solution to the equation Ax = b into ans1, ans2, ans3, where A is 3x3 and b is 3x1. (Function solves system using cramers rule. )
 
 minimumIndex - Returns the index corresponding to the minimum of an array.
 
 copyFloatArrayToDoubleArray - Converts a float array to double array. 
 
 
 
Definition of parameters: 

 double wavelengths[] - array of wavelengths. (Independent variable)
 
 double data[] - spectrum data corresponding to wavelengths
 
 int numberOfChannels - number of channels in spectrum
 
 int lowerIndex - channel number of band lower boundary 
 
 int upperIndex - channel number of band upper boundary 
 
 int minIndex - channel number of band minimum 
 
 double bandmin - wavelength of band minimum
 
 double bandDepth - difference between continuum value and minimum band value
 
 double FWHMx_,FWHMy_ - the interpolated coordinate of half max point
 
 double c2,c1,c0 - coefficients of fitted parabola in the form: c2x^2 + c1x + c0
 
 
 In all cases: lowerIndex < minIndex < upperIndex < length of wavelengths & data
------------------------------------------------------------------------------*/ 

// ******************* Function Prototypes *******************

// function: calcAndDisplayBandParameters2
// Given a spectrum and a list of continuum boundaries as ranges, the band parameters are calculated and printed to the command line for each band. 
void calcAndDisplayBandParameters2(double wavelengths[], double data[], double lowerbound1[], double upperbound1[], double lowerbound2[], double upperbound2[], int numberOfChannels, int numberOfBands);

// function: calcAndDisplayBandParameters
// Given a spectrum and a list of continuum boundaries as points, the band parameters are calculated and printed to the command line for each band. 
void calcAndDisplayBandParameters(double wavelengths[], double data[], double lowerbound[], double upperbound[], int numberOfChannels, int numberOfBands);

// function: calcParameters2
// Saves band minimum, band depth, FWHM, and Asymmetry to pointers passed as parameters. d1 and d2 are the continuum data values corresponding to the lowerIndex and upperIndex. Used for continuum boundaries passed as ranges. i.e. the continuum value is not necessarily a measured value.  
void calcParameters2(double wavelengths[], double data[], double d1, double d2, int lowerIndex, int upperIndex, int numberOfChannels, double *bandMin, double *bandDepth, double *FWHM, double *asym);

// function: calcParameters
// Saves band minimum, band depth, FWHM, and Asymmetry to pointers passed as parameters. 
void calcParameters(double wavelengths[], double data[], int lowerIndex, int upperIndex, int numberOfChannels, double *bandMin, double *bandDepth, double *FWHM, double *asym);

// Given an array of wavelengths, corresponding data, and arrays of continuum point boundaries, calculates the continuum points for those ranges. 
// upperbool = 1 if the range is around the upper continuum pt. 
// Returns the index of the channel and the continuum value.   
void findCpoint(double wavelengths[], double data[], double bounds1[], double bounds2[], int ci[], double cval[], int numberOfBands, int numberOfChannels, int upperbool);

// Produces the elements of channels[] which correspond to the nearest channels to the contiuum points in bounds[] 
void matchToChannels(double A[], double bounds[], int channels[], int numberOfBands, int numberOfChannels );

// Searchs through array, A, and finds the element nearest to the target value. 
// Returns the index of that value.
int channelSearch( double A[], double target, int F, int L);

// function: calcAsymmetry
// Returns the asymmetry of the specified band.
double calcAsymmetry(double wavelengths[], double data[], int lowerIndex, int upperIndex, int minIndex, double bandDepth , double bandmin );

// function: removeContinuum
// Saves the continuum removed data to newdata[] (which is passed by user.) 
void removeContinuum(double wavelengths[], double data[], double newdata[], double d1, double d2, int lowerIndex, int upperIndex, int numberOfChannels);

// function: calcFWHM
// Returns the Full Width at Half Maximum of band. 
double calcFWHM (double wavelengths[], double data[], int lowerIndex, int upperIndex, int minIndex, double bandDepth);

// function: calcFWHMpts
// Saves the two coordinates interpolated to fall on the half max value. 
void calcFWHMpts( double wavelengths[], double data[], int lowerIndex, int upperIndex, int minIndex, double bandDepth, double *FWHMx1, double *FWHMy1, double *FWHMx2, double *FWHMy2);

// function: calcBandDepth 
// Returns the difference between continuum value and minimum band value.
double calcBandDepth(double wavelengths[], double data[], int lowerIndex, int upperIndex);

// function: calcBandMin
// Returns the wavelength at which the band is minimum based on parabolic fit. 
double calcBandMin(double wavelengths[], double data[], int lowerIndex, int upperIndex );

// function: getParabolicFit
// Fits a parabola to the lowest three data points and saves the coefficients in coeff2, coeff1, coeff0.
void getParabolicFit(double wavelengths[], double data[], int minIndex, double *coeff2, double *coeff1, double *coeff0);

// function: parabolaEval
// Evaluates a parabola function of the form c2x^2 + c1x + c0 @ xval
double parabolaEval( double c2, double c1, double c0, double xval );

// function: calcparabolaMinx
// Returns x-coordinate of minimum of parabolic function 
// c2, c1, c0 represent a parabola of the form c2x^2 + c1x + c0
double calcparabolaMinx( double c2, double c1, double c0 );

// function: determinant3by3
// Returns the determinant of 3x3 matrix
double determinant3by3( double A[3][3] );

// function: copy3by3
// Copies 3x3 matrix "in" to 3x3 matrix "out"
void copy3by3( double in[3][3], double out[3][3]);

// function: solve3by3system
// Assigns the solution to the equation Ax = b into ans1, ans2, ans3
// where A is 3x3 and b is 3x1. (Function solves system using cramers rule. )
void solve3by3system(double A[3][3], double b[3], double *ans1, double *ans2, double *ans3);

// function: minimumIndex
// Returns the index corresponding to the minimum of an array.
int minimumIndex(double array[], int lowerindex, int upperindex);

// function: copyFloatArrayToDoubleArray
// Converts and array of float values to double values. 
//void copyFloatArrayToDoubleArray( float input[] , double output[], int size);



// ******************* Function Definitions *******************

void calcAndDisplayBandParameters2(double wavelengths[], double data[], double lowerbound1[], double upperbound1[], double lowerbound2[], double upperbound2[], int numberOfChannels, int numberOfBands)
{
	// indicies of channels corresponding to the continuum points.
	int lowerCpointi[numberOfBands], upperCpointi[numberOfBands];
	
	// y-values of the continuum points. (Average of values within range)
	double lowerCpointval[numberOfBands], upperCpointval[numberOfBands];
	
	// Declare placeholders for band parameters of interest. 
	double bandMin; 	// Minimum data point in the band range. 
	double bandDepth; 	// Depth of the band.
	double FWHM; 		// The full width at half maximum of the band. 
	double asym; 		// The calculated asymmetry of the band. 
	int i;				// loop index
	
	// Fills the lowerCpointi and lowerCpointval arrays
	findCpoint(wavelengths, data, lowerbound1, upperbound1, lowerCpointi, lowerCpointval, numberOfBands, numberOfChannels, 0);
	
	// Fills the upperCpointi and upperCpointval arrays
	findCpoint(wavelengths, data, lowerbound2, upperbound2, upperCpointi, upperCpointval, numberOfBands, numberOfChannels, 1);

	// Calculates band min, band depth, fwhm, asymmetry.
	for(i = 0; i < numberOfBands; i++)
	{
		printf("Band Wavelengths: %f - %f\n",wavelengths[lowerCpointi[i]],wavelengths[upperCpointi[i]]);
		printf("Band Channels: %d - %d\n",lowerCpointi[i] + 1,upperCpointi[i] + 1);
		
		calcParameters2(wavelengths, data, lowerCpointval[i], upperCpointval[i], lowerCpointi[i], upperCpointi[i], numberOfChannels, &bandMin, &bandDepth, &FWHM, &asym);
	
		// Print results to command line. 
		printf("Band Min\tBand Depth\tF.W.H.M.\tAsymmetry\n");
		printf("%f\t%f\t%f\t%f\n\n",bandMin,bandDepth,FWHM,asym);
		
	}
}

void calcAndDisplayBandParameters(double wavelengths[], double data[], double lowerbound[], double upperbound[], int numberOfChannels, int numberOfBands)
{
	int lowerindex[numberOfBands], upperindex[numberOfBands];
	// Declare placeholders for parameters of interest. 
	double bandMin; 	// Minimum data point in the band range. 
	double bandDepth; 	// Depth of the band.
	double FWHM; 		// The full width at half maximum of the band. 
	double asym; 		// The calculated asymmetry of the band. 
	int i;
	
	// Convert input arrays of wavelength boundary to arrays of channel numbers.
	matchToChannels( wavelengths, lowerbound, lowerindex, numberOfBands, numberOfChannels);
	matchToChannels( wavelengths, upperbound, upperindex, numberOfBands, numberOfChannels);
	
	// Calculates band min, band depth, fwhm, asymmetry.
	for(i = 0; i < numberOfBands; i++)
	{
		printf("Band Wavelengths: %f - %f\n",wavelengths[lowerindex[i]],wavelengths[upperindex[i]]);
		printf("Band Channels: %d - %d\n",lowerindex[i] + 1,upperindex[i] + 1);
		
		calcParameters(wavelengths, data, lowerindex[i], upperindex[i], numberOfChannels, &bandMin, &bandDepth, &FWHM, &asym);
		
		// Print results to command line.
		printf("Band Min\tBand Depth\tF.W.H.M.\tAsymmetry\n");
		printf("%f\t%f\t%f\t%f\n\n",bandMin,bandDepth,FWHM,asym);
		 
	}
}



void calcParameters2(double wavelengths[], double data[], double d1, double d2, int lowerIndex, int upperIndex, int numberOfChannels, double *bandMin, double *bandDepth, double *FWHM, double *asym)
{
	double data2[numberOfChannels];	// The continuum removed data. 
	int minIndex; 		// The index corresponding to the band minimum.
	// bandMin - Minimum data point in the band range. 
	// bandDepth - Depth of the band.
	// FWHM - The full width at half maximum of the band. 
	// asym - The calculated asymmetry of the band. 
	
	// Remove the continuum from band
	removeContinuum(wavelengths, data, data2, d1, d2, lowerIndex, upperIndex, numberOfChannels);
	
	// Search for index of band minimum.
	minIndex = minimumIndex(data2, lowerIndex, upperIndex);
	
	// Calculate the band minimum based on parabolic fit of lowest three data points. 
	*bandMin = calcBandMin(wavelengths, data2, lowerIndex, upperIndex );	
	
	// Calculate the band depth by differencing the continuum and band minimum.
	*bandDepth = calcBandDepth(wavelengths,data2, lowerIndex, upperIndex);
	
	// Calculate the full width at half maximum.
	*FWHM = calcFWHM(wavelengths, data2, lowerIndex, upperIndex, minIndex, *bandDepth );
	
	// Calculate the asymmetry of the band.
	*asym = calcAsymmetry(wavelengths, data2, lowerIndex, upperIndex, minIndex, *bandDepth, *bandMin);

}


void calcParameters(double wavelengths[], double data[], int lowerIndex, int upperIndex, int numberOfChannels, double *bandMin, double *bandDepth, double *FWHM, double *asym)
{
	double data2[numberOfChannels];	// The continuum removed data. 
	int minIndex; 		// The index corresponding to the band minimum.
	// bandMin - Minimum data point in the band range. 
	// bandDepth - Depth of the band.
	// FWHM - The full width at half maximum of the band. 
	// asym - The calculated asymmetry of the band. 
	
	// Remove the continuum from band
	removeContinuum(wavelengths, data, data2, data[lowerIndex], data[upperIndex] , lowerIndex, upperIndex, numberOfChannels);
	
	// Search for index of band minimum.
	minIndex = minimumIndex(data2, lowerIndex, upperIndex);
	
	// Calculate the band minimum based on parabolic fit of lowest three data points. 
	*bandMin = calcBandMin(wavelengths, data2, lowerIndex, upperIndex );	
	
	// Calculate the band depth by differencing the continuum and band minimum.
	*bandDepth = calcBandDepth(wavelengths,data2, lowerIndex, upperIndex);
	
	// Calculate the full width at half maximum.
	*FWHM = calcFWHM(wavelengths, data2, lowerIndex, upperIndex, minIndex, *bandDepth );
	
	// Calculate the asymmetry of the band.
	*asym = calcAsymmetry(wavelengths, data2, lowerIndex, upperIndex, minIndex, *bandDepth, *bandMin);

}

// Given an array of wavelengths, corresponding data, and arrays of continuum point boundaries, calculates the continuum points for those ranges. 
// upperbool = 1 if the range is around the upper continuum pt. 
// Returns the index of the channel and the continuum value.   
void findCpoint(double wavelengths[], double data[], double bounds1[], double bounds2[], int ci[], double cval[], int numberOfBands, int numberOfChannels, int upperbool)
{
	// Holds the channel index corresponding to the boundary wavelengths.  
	int bounds1i[numberOfBands], bounds2i[numberOfBands];
	// i,j are indicies and count is the count of data points to divide by when finding the average data point between the bounds.  
	int i, j, count;
	
	// Fills the bounds1i array with channel indicies of bounds1
	matchToChannels( wavelengths, bounds1, bounds1i, numberOfBands, numberOfChannels );
	// Fills the bounds2i array with channel indicies of bounds2
	matchToChannels( wavelengths, bounds2, bounds2i, numberOfBands, numberOfChannels );
	
	// Iterate through all of the bands and find the channel index and value of continuum point. 
	for( i = 0; i < numberOfBands; i++ )
	{
		count = 0;
		// The continuum index is the average between the boundary indicies.
		ci[i] = (bounds1i[i] + bounds2i[i]) / 2;
		// If the average is in between two channels, take the outer one.
		if( bounds1i[i] + bounds2i[i] % 2 != 0 || upperbool == 1)
		{
			ci[i] += 1;
		}
			
		// Find the average y-val and assign it to cval. 	
		for( j = bounds1i[i]; j <= bounds2i[i]; j++ )
		{
			cval[i] += data[j];
			count++;
		}
		cval[i] /= count;
	}
	
}

// Given a list of wavelengths 
// Returns the indicies of the "closest" channel in the input spectrum.
void matchToChannels(double A[], double bounds[], int channels[], int numberOfBands, int numberOfChannels )
{
	double target;
	int i;
	for( i = 0; i < numberOfBands; i++ )
	{
		target = bounds[i];
		channels[i] = channelSearch(A,target,0,numberOfChannels - 1);
	}
}

// Searchs through array, A, and finds the element nearest to the target value. 
// Returns the index of that value. 
int channelSearch( double A[], double target, int F, int L)
{
	int i, channel = 0;
	double diff, mindiff;
	mindiff = fabs(A[F] - target);
	
	for( i = F + 1; i <= L; i++)
	{
		diff = fabs(A[i] - target);
		if(diff < mindiff)
		{
			mindiff = diff;
			channel = i;
		}
	}
	return channel;
}

double calcAsymmetry(double wavelengths[], double data[], int lowerIndex, int upperIndex, int minIndex, double bandDepth , double bandmin )
{
	double FWHMx1,FWHMy1,FWHMx2,FWHMy2;
	calcFWHMpts(wavelengths, data, lowerIndex, upperIndex, minIndex, bandDepth, &FWHMx1,  &FWHMy1,  &FWHMx2, &FWHMy2);
	return fabs( (FWHMx1 - bandmin) / (FWHMx2 - bandmin) );
}

void removeContinuum(double wavelengths[], double data[], double newdata[], double d1, double d2, int lowerIndex, int upperIndex, int numberOfChannels)
{
	double w1,w2,slope;
	int i;
	w1 = wavelengths[lowerIndex];
	w2 = wavelengths[upperIndex];
	slope = (d2 - d1)/(w2 - w1);
	for(i = 0; i < numberOfChannels; i++)
		newdata[i] = data[i]/(slope*(wavelengths[i] - w1) + d1);
}

double calcFWHM (double wavelengths[], double data[], int lowerIndex, int upperIndex, int minIndex, double bandDepth)
{
	double FWHMx1,FWHMy1,FWHMx2,FWHMy2;
	calcFWHMpts(wavelengths, data, lowerIndex, upperIndex, minIndex, bandDepth, &FWHMx1,  &FWHMy1,  &FWHMx2, &FWHMy2);
	return FWHMx2 - FWHMx1;
}

void calcFWHMpts( double wavelengths[], double data[], int lowerIndex, int upperIndex, int minIndex, double bandDepth, double *FWHMx1, double *FWHMy1, double *FWHMx2, double *FWHMy2)
{
	double halfmax;
	double x1,x2,y1,y2;
	
	
	// Pull out the region of interest
	double diffsleft[minIndex - lowerIndex], diffsright[upperIndex - minIndex];
	
	int i,leftIndex1,leftIndex2,rightIndex1,rightIndex2;
	
	
	// Find the half max value.
	halfmax = 1 - bandDepth/2.0;
		
	//Find the FWHM points based on a linear interpolation between the two nearest data points in the spectrum.

	//Search for the closest intensity to the half max on both sides of the band. 
	//Record the index at which this occurs.
	
	for( i = 0; i < minIndex - lowerIndex; i++ )
	{
		diffsleft[i] = fabs(data[lowerIndex + i] - halfmax);
	} 
	leftIndex1 = minimumIndex(diffsleft,0,minIndex - lowerIndex - 1) + lowerIndex;
	
	// Find the next closest to linearly interpolate over the region of interest. 
	if( fabs(data[leftIndex1 - 1] - halfmax) < fabs(data[leftIndex1 + 1] - halfmax))
	{
		leftIndex2 = leftIndex1;
		leftIndex1--; 
	}
	else
		leftIndex2 = leftIndex1 + 1;
		
	// Perform the linear interpolation to achieve the appropriate FWHM point
	x1 = wavelengths[leftIndex1];
	y1 = data[leftIndex1];
	x2 = wavelengths[leftIndex2];
	y2 = data[leftIndex2];
	*FWHMx1 = (halfmax - y1)*(x2 - x1)/(y2-y1) + x1;
	*FWHMy1 = (y2-y1)/(x2-x1)*(*FWHMx1 - x1) + y1;

	// And now for the right side...
	for( i = 0; i < upperIndex - minIndex; i++ )
	{
		diffsright[i] = fabs(data[minIndex + i] - halfmax);
	} 
	rightIndex1 = minimumIndex(diffsright,0,upperIndex - minIndex - 1);
	rightIndex1 += minIndex;
	
	// Find the next closest to linearly interpolate over the region of interest. 
	if( fabs(data[rightIndex1 - 1] - halfmax) < fabs(data[rightIndex1 + 1] - halfmax))
	{
		rightIndex2 = rightIndex1;
		rightIndex1--; 
	}
	else
		rightIndex2 = rightIndex1 + 1;
		
	// Perform the linear interpolation to achieve the appropriate FWHM point
	x1 = wavelengths[rightIndex1];
	y1 = data[rightIndex1];
	x2 = wavelengths[rightIndex2];
	y2 = data[rightIndex2];
	*FWHMx2 = (halfmax - y1)*(x2 - x1)/(y2-y1) + x1;
	*FWHMy2 = (y2-y1)/(x2-x1)*(*FWHMx2 - x1) + y1;
}

// Assumes the continuum value is 1.
double calcBandDepth(double wavelengths[], double data[], int lowerIndex, int upperIndex)
{
	int minIndex;
	double c2,c1,c0,bandMinY; // parabola of form: c2x^2 + c1x + c0
	minIndex = minimumIndex(data,lowerIndex,upperIndex);
	getParabolicFit(wavelengths,data,minIndex,&c2,&c1,&c0);
	bandMinY = parabolaEval(c2,c1,c0,calcparabolaMinx(c2,c1,c0));
	return 1 - bandMinY;
}

double calcBandMin(double wavelengths[], double data[], int lowerIndex, int upperIndex )
{
	int minIndex;
	double c2,c1,c0; // parabola of form: c2x^2 + c1x + c0
	minIndex = minimumIndex(data,lowerIndex,upperIndex);
	getParabolicFit(wavelengths,data,minIndex,&c2,&c1,&c0);
	return calcparabolaMinx(c2,c1,c0);
}

void getParabolicFit(double wavelengths[], double data[], int minIndex, double *coeff2, double *coeff1, double *coeff0)
{
	double A[3][3];
	double b[3];
	A[0][0] = wavelengths[minIndex - 1]*wavelengths[minIndex - 1];
	A[1][0] = wavelengths[minIndex]*wavelengths[minIndex];
	A[2][0] = wavelengths[minIndex + 1]*wavelengths[minIndex + 1];
	A[0][1] = wavelengths[minIndex - 1];
	A[1][1] = wavelengths[minIndex];
	A[2][1] = wavelengths[minIndex + 1];
	A[0][2] = 1;
	A[1][2] = 1;
	A[2][2] = 1;
	b[0] = data[minIndex - 1];
	b[1] = data[minIndex];
	b[2] = data[minIndex + 1];
	solve3by3system(A,b,coeff2,coeff1,coeff0);
}

double parabolaEval( double c2, double c1, double c0, double xval )
{
	return c2*xval*xval + c1*xval + c0;
}

double calcparabolaMinx( double c2, double c1, double c0 )
{
	return (-1.0*c1)/(2.0*c2);
}

// Returns the determinant of 3x3 matrix
double determinant3by3( double A[3][3] )
{
	double det = 0;
	det += A[0][0]*(A[1][1]*A[2][2] - A[1][2]*A[2][1]);
	det -= A[0][1]*(A[1][0]*A[2][2] - A[1][2]*A[2][0]);
	det += A[0][2]*(A[1][0]*A[2][1] - A[1][1]*A[2][0]);
	return det;
}

// copies 3x3 matrix "in" to 3x3 matrix "out"
void copy3by3( double in[3][3], double out[3][3])
{
	int i,j;
	for( i = 0; i < 3; i++)
		for( j = 0; j < 3; j++)
			out[i][j] = in[i][j];
}


// assigns the solution to the equation Ax = b into ans1, ans2, ans3
// where A is 3x3 and b is 3x1. Solves system using cramers rule. 
// Example usage:
//	double c2,c1,c0; // form: c2x^2 + c1x + c0
//	solve3by3system(a,b,&c2,&c1,&c0);
//	printf("%f\n%f\n%f\n",c2,c1,c0); 
void solve3by3system( double A[3][3], double b[3], double *ans1, double *ans2, double *ans3 )
{
	double D, D1, D2, D3;
	double A1[3][3], A2[3][3], A3[3][3];
	copy3by3(A,A1);
	copy3by3(A,A2);
	copy3by3(A,A3);
	A1[0][0] = b[0];
	A1[1][0] = b[1];
	A1[2][0] = b[2];
	A2[0][1] = b[0];
	A2[1][1] = b[1];
	A2[2][1] = b[2];
	A3[0][2] = b[0];
	A3[1][2] = b[1];
	A3[2][2] = b[2];
	D = determinant3by3(A);
	D1 = determinant3by3(A1);
	D2 = determinant3by3(A2);
	D3 = determinant3by3(A3);
	*ans1 = D1/D;
	*ans2 = D2/D;
	*ans3 = D3/D;
}

// Returns the the index corresponding to the minimum of an array
int minimumIndex(double array[], int lowerindex, int upperindex)
{
	int i,iMin;
	double min;
	min = array[lowerindex];
	iMin = lowerindex;
	for( i = lowerindex + 1; i <= upperindex; i++ )
	{
		if( array[i] < min )
		{
			iMin = i;
			min = array[i];
		}
	}
	return iMin;
}

//void copyFloatArrayToDoubleArray( float input[] , double output[], int size)
//{
//	int i;
//	for(i = 0; i < size; i++)
//		output[i] = input[i];
//}
